c
      subroutine build_sx(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
      real*8 Ip(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o,Ip
      real*8 ms
      real*8 lm(n)

      Ip = 0
      o=0

      lm=0
      ms=-(n-1)/2.0
      do k=1,n,1
      lm(k) = ms
      ms=ms+1
      end do

      do k=1,n-1,1
      Ip(k , k + 1) = sqrt(((n-1)/2.0 - lm(k)) * ((n-1)/2.0 + lm(k)+ 1))
      end do

      o = (1.0 / 2) * (Ip + transpose(Ip))
c      return
      end subroutine build_sx
c
c
c
      subroutine build_sy(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
      real*8 Ip(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o,Ip
      real*8 ms
      real*8 lm(n)

      Ip = 0
      o=0

      lm=0
      ms=-(n-1)/2.0
      do k=1,n,1
      lm(k) = ms
      ms=ms+1
      end do

      do k=1,n-1,1
      Ip(k , k + 1) = sqrt(((n-1)/2.0 - lm(k)) * ((n-1)/2.0 + lm(k)+ 1))
      end do

      o = (1.0 / 2) * (-Ip + transpose(Ip))
c      return
      end subroutine build_sy
c
c
      subroutine build_sz(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o
      real*8 ms

      o=0

      ms=(n-1)/2.0
      do k=1,n,1
      o(k,k)=ms
      ms=ms-1
      end do
c      return
      end subroutine build_sz
c
c
      subroutine build_sd(o,n)
      implicit none
c a subroutine to build an identity matrix with a dimension n
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o
      o=0
      do k=1,n,1
      o(k,k)=1
      end do
c      return
      end subroutine build_sd
c
c
      subroutine bfield(o,b0,xyz,f,phi,A,t)
      real*8 o(3)
Cf2py intent(out) o
      real*8 xyz(3)
Cf2py intent(in) xyz
      real*8 b0,f,phi,A,t
Cf2py intent(in) b0,f,phi,A,t
      real*8 pi,norm

      pi=3.1415026535
      o=0
      norm=sqrt(xyz(1)**2+xyz(2)**2+xyz(3)**2)
      o = (b0 + A* sin(2*pi*t*f+phi)) / norm * xyz
c      return
      end subroutine bfield
c
c


